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Abstract 

Two variants of the MCV3 scheme are presented based on a flux reconstruction formulation. Different from the 
original multi-moment constrained finite volume method of third order (MCV3), the multi-moment constraints are 
imposed at the cell center on the point value, the first and second order derivatives. The continuity of the flux function 
at cell interfaces are also used as the constraints to ensure the numerical conservation. Compared to the original 
MCV3 scheme, both two variants have higher numerical accuracy and less restrictive CFL condition for computational 
stability. Moreover, without the need to solve derivative Riemann problem at cell boundaries, the new schemes benefit 
the implementations in arbitrary quadrilateral in 2D and hexahedron in 3D. 
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1. The MCV3 scheme 

We consider the following conservation law 

Q-f where u is the solution function, and f(u) the flux function. The computational domain is divided into / non- 
overlapping cells or elements Q, = [x^i, x i+ i], i = 1,2,- •• ,1. IntheMCV3 scheme 3 solution points x±,k- 1,2,3, 
are set over £2, where the solution u$-, k = 1,2,3, is computed. Suppose that a proper approximation of the flux 
function fi(x) is reconstructed, we can immediately update the solutions within £2, by the following point-wise semi- 
discretized equations at solution points x^, 



dt 



dfi(x) 



dx 



k= 1,2,3. (2) 



The central task left now is how to reconstruct the flux function fi(x). In principle, the way to reconstruct /j(jc) 
makes difference among the numerical schemes|[ll[3t]. 

For brevity, we make use of a local coordinate system £ e [—1,1] that transforms the real mesh cell x e 

[Xi-i/2,Xi+i/ 2 ] by 



X - Xji 



€ = 2—^.-1, (3) 

Axt 



where Ajc,- = Xj+i/2 - Xi-\/2 is the mesh spacing. 

The time evolution equations for updating the solutions are correspondingly 
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We assume that the flux f(u) is a function of solution u. The values of the flux at the corresponding points, fa, 
k- 1, 2, 3, are computed directly. The primary reconstruction of flux function is then built by 



fm = Y J fa<Pk{&. 



k=l 



where 



n 



(5) 



(6) 



is the basis function of the Lagrange interpolation. Interpolation (f5J is constructed piecewisely over each cell element 
without connection to its neighboring cells, thus cannot be directly used to update the solutions. A modification is 
required. 

In the MCV3 scheme, the modified flux function is constructed by imposing the continuities of the flux function 
and its first order derivative at the cell boundaries, which is realized by the following constraint conditions, 



where 



C/k-i) = if(-i); 
Mi) = /f (i); 
4/ ] (-i) = 4/ ls (-i); 

4/ ] (i) = /£ ]S (d, 



//(-D = /f (-1) = /*,, 
Ml) = f(D = Jti 



(7) 



(8) 



and 



4- ] (-i)=4r(- 1 ) 



f[l]S/ 



[1]S 



(9) 



1|S 



The numerical flux at the cell boundary x t _ i , , , and its first order derivative, are solved by following Riemann 
problems(DRP), 



ft 



Riemann 



/^^DRtemann^,^), 



(10) 



where "Riemann(- , ■)" and "DRiemann(- , •)" denote the solvers for the conventional and derivative Riemann prob- 
lems respectively. 

(|7]i is a Hermite interpolation to determine the modified flux function which is written in a polynomial form as, 



+\(f§m-f§(-v)e 

+ l - (3/f (i) - 3/f(-i) - /|(-i) - 

+ T (2/i S (D + 2f'f{-\) + - /|(1)) . 



(11) 



The first order derivative (gradient) of (0 reads then, 

3 



fdo =4 (/f (-d - fm + /|(-d + /|(d)^ 2 

+5(jf(D-Jf(-D)f 
4( 3 ^ 1 )- 3 ^- 1 )-/^- 1 )-/^ 1 ))- 



(12) 



The derivatives of the modified flux function at the solution points are obtained as 



(13) 



v I df A3 

The solutions are then immediately computed by © with a proper time integration algorithm. 

In the original MCV3 scheme[2], the solution points are equally spaced and including two cell ends, i.e. £i = — 1, 
^2=0 and £3 = 1. The left/right-most solution points coincide with the cell boundaries. In this case, the continuity 
conditions of flux function at the cell boundaries are automatically satisfied, and only the derivatives of the flux 
function need to be computed from the DRP. 

The derivatives of the modified flux function at the solution points are obtained as 

\fm =/ f i«i) = /!(-i); 

hi = \ (3./f(D - 3/f (-1) - - /|(1)) ; (14) 

l/fi3=/fi(ft)=/|(l). 

It is straightforward to show the following conservation property, 



k=l v 



1 r 



4 - 



1 ^ 



(15) 



The solution points can be chosen as other quadrature point sets, such as the Legendre or Chebyshev Gauss points, 
but we find from Fourier analysis and numerical tests that the different solution point sets don't alter significantly the 
numerical result. 



2. The variants of MCV3 scheme 



We present here two variants by making use of different constraints in determining the modified flux function. 
Instead of the constraint conditions of (0, we impose the multi-moment constraints at the cell center, 

f/k-D = /f(-D; 

Md = fay, 

fi(0) = MO); (16) 

/]/ 1] (°) = 4i ] (°); 
./? ] (0) = 4 2I (0)- 

We retain the continuity of the modified flux function at cell boundaries, which is required from the numerical con- 
servation and stability. The rest of the constraints are determined from the primary interpolation function in terms of 
derivatives. 
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2.1. MCV3 scheme for uniform points with center constraints: MCV3JJPCC 

Same as in the original MCV3 scheme[2], the solution points are equally spaced and including two cell ends, 
ft = -1,6=0 andft = 1. 

Constraint conditions ([Tol l allows to reconstruct a polynomial of 4th degree, 



7;(ft =~ (/f (-D + /f(D - fa - fa)? + ~ (/ S (D - /f(-D + /i - fa)? 



The first-order derivative then reads, 



(17) 



/fi(ft =2 (/f (-1) + ,/f(l) - / fl - h)? + 2 (/f d) - ./f (-1) + /n - fo)f 

+ (fn-2fi 2 +fi3){+\(fi3-fi). 
The derivatives of the modified flux function at the solution points are obtained as 

1 — „<R , - * 



/f'l 


= /fKft) 




= /f,(6) 


fen 


= /?;(ft) 



rC/b - /ii); 

„ . 1 



It is straightforward to show the following conservation property, 

2 [fm £ = \h + \ha + \fo = ff{\) - ff{-D- 



(18) 



(19) 



(20) 



2.2. MCV3 scheme for Chebyshev points with center constraints: MCV3J2PCC 

We use the Chebyshev-Gauss points, i.e. ft = - v3/2, ft = and ft = v3/2, as the solution points. Constraint 
conditions (fTol i leads to the following polynomial of 4th degree, 



/■(ft =i (3/f(-D + 3/f (1) - 4/i, + 2/ Q - 4/e)* 4 

l - [iff {I) - 3/f (-1) + 2 V§/ a - 2 V3/b)£ 3 
V3 



+t (4/h - 8/ a + 4/ i3 ) ft + — (/; 3 - A) f + / a . 
O 3 



The first-order derivative then reads, 



(21) 



//(ft =3 (3/ S (-D + 3/f (1) - A f a + 2 f i2 - 4f i3 )f 

+ \ (3/ S (D - 3/ s (-l) + 2 V3/ a - 2 V3/i 3 )ft 

1 V3 
+ 3 (4/h - 8/ a + 4/ s ) ? + — (fc - /a) • 

The derivatives of the modified flux function at the solution points are obtained as 



fm 


= /f,(ft) 


ha. 


= /fi(ft) 


fm 


= /f-(ft) 



3 /3 



/fa = /f«(ft) = j V3/i + z V3/a - ^ V3/ i3 + ~ ( - - V3~)/f (1) - ± ( ± + V3)^(-l); 



412 



3 (3 
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V3 



(/■3-/i); 



/fa = /f,(ft) = ^ V3/ a - I V3/a - ~ V3/ i3 + + V3~)/f (1) Vaj/f (-1). 



(22) 



(23) 



From ( 1231 . the numerical conservation can be immediately proved by the following equality, 

2 [f&k £ = pw + + ^/ f /3 = /fd) - f'f(-l)- (24) 



Semidiscrete Spectra Semidiscrete Spectra Semidiscrete Spectra 




Figure 1: The spectrum of the semi-discrete schemes. Left: MCV3; Center: MCV3.UPCC, Right: MCV3.CPCC. 



3. Fourier analysis 

In this section, we evaluated the numerical schemes previously discussed by examining the Fourier mode trans- 
ported with the following advection equation. 

g + ^=0. (25) 
at ox 

Using a wave solution 

q(x,t) = e I( * x * t \ (26) 

and assuming a uniform grid spacing x i+ i - x t _i = Ax, we have = e /w O,+f/A-v/2) an( j M(; = e _/wA *H ijt , which recast 
the time evolution equations for the solutions into 

§=Su, (27) 
at 

The properties of the numerical schemes can be examined by analyzing the eigenvalues of ( 1271 ). FigQ]shows the 
spectrum (collection of all eigen values) of S for different schemes. It is observed that all eigenvalues lie on the left 
half of the real axis, i.e the negative real parts indicate that all the schemes are stable under the CFL conditions. The 
allowable CFL numbers for computational stability can be estimated by the largest eigenvalue, the spectral radius 
p for each scheme, i.e. a scheme has a larger spectral radius has to use a smaller CFL number for computational 
stability. We know from Fourier analysis that Pmcv3.upcc = 4.7, Pmcv3.cpcc =5.4 and Pmcv3 = 6.0, which reveals 
that MCV3_CPCC scheme has the largest stable CFL number. This is confirmed by numerical tests for the linear 
advection equation. With a 3rd-order Runge Kutta scheme, the largest allowable CFL numbers are 0.47, 0.44 and 
0.41 forMCV3.UPCC, MCV3_CPCC and MCV3 respectively. 
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Numerical dispersion 




Figure 2: Numerical dispersion relations for four-point schemes. 



The numerical errors of different schemes can be examined by comparing the principal eigenvalue of S, ^g (LllA( y 
with the exact solution, -Ioj, of the advection equation d25l l for initial condition, 

q(x,0) = e Iwx . (28) 



The error of a given semi-discrete formulation is 



and the convergence rate is evaluated by 



w= K^I) /in(2) )" L (3o) 



Table 1: Numerical errors and convergence rates. 



Scheme 


ojAx 


= 7T/8 


toAx 


= tt/16 


order 


MCV3 


-3.25 x 1(T 4 


-3.33x 10-=/ 


-2.06 x KT 5 


- 1.07 x 10- b i 


2.99 


MCV3.UPCC 


-1.65 x 1(T 4 


- 2.22 x 1(T 6 / 


-1.03 x 10~ 5 


-6.81 x 10~ 9 / 


3.00 


MCV3.CPCC 


-3.15 x 1(T 6 


+ 1.91 x !Q- 5 i 


-4.93 x 10~ 8 


+ 6.05 x 10~ 7 ; 


3.99 



The numerical errors of the three schemes are given in Table 1 . MCV3_UPCC and MCV3_CPCC are more accurate 
than the original MCV3 scheme. Similar to the original MCV3, MCV3JUPCC shows a 3rd-order convergence rate, 
while MCV3_CPCC has a 4th-order convergence rate. 

The dispersion and dissipation relations of the spatial discretization can be evaluated by plotting the real and 
imaginary parts of the principal eigenvalues as a functions of the wave number a>. From Figj2] we find that the 
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Numerical dissipation 
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Figure 3: Numerical dispersion relations for three-point schemes. 



MCV3JJPCC has the most accurate numerical dispersion. MCV3_CPCC has a different dispersion behavior com- 
pared to MCV3 and MCV3JJPCC. The numerical dissipations are plotted in FigfJ] MCV3_CPCC is the best, which 
is also observed in numerical tests. Both MCV3_UPCC and MCV3_CPCC have improved dissipation accuracy over 
the original MCV3 scheme. 

We further give the Taylor expansion of the eigen values of S with respect to the mesh size in Table 2. 



Table 2: Taylor expansion of the eigenvalues in terms of the mesh size. The principal eigenvalue ^g^/^ is Ai for all schemes. 



Scheme 


Bit 


ien values 


MCV3 


M 


= -Ico - 1.39 x 10^w 4 Ajc 3 - 3.70 x lO^I^Ax 4 + O(A^) 




h 


= -6/ Ax + 3ICJ + to 2 Ax 




A 3 


= 


MCV3.UPCC 


M 


= -1(0 - 6.94 x KTWAji 3 - 2.31 x lO^/VAx 4 + OiAx?) 

= (-4 + 2.83/)/ A* + (0.71 + T)(0 + (0.25 - 8.84 2 /)w 2 Ax + 0(Ax 2 ) 




A 2 




A 3 


= (-4 - 2.83/)/Ax + (-0.71 + I)oj + (0.25 + 8.84 2 /)w 2 Ax + 0(Ax 2 ) 


MCV3_CPCC 




= -Ico + 2.08 x lO^/VA* 4 - 8.68 x 10- 4 w b Ax 5 + 0(Ax b ) 




h 


= (-3 + 3.87/)/Ax + (1.16 + 0.5I)oj - 0.31/w 2 Ajc + 0(A^ 2 ) 




A 3 


= (-3 - 3.87/)/Ax + (-1.16 + 0.5/)w + 0.31/w 2 Ax + 6>(Ax 2 ) 



Consistent with the observations aforementioned, we find from Table 2 that MCV3JUPCC and MCV3_CPCC 
have less truncation errors in both dissipation (real part) and dispersion (imaginary part) compared to the original 
MCV3. MCV3 and MCV3_UPCC have a third order accuracy in dissipation and all schemes have a fourth order 
accuracy in dispersion. The MCV3_CPCC is superior in dissipation accuracy which is fifth order, two orders higher 
than the others. The computational modes are represented by A 2 and A 3 . It is observed that all the schemes have real 
negative parts with the leading terms of order 0(Ax~ l ) in A 2 and A 3 , which means that the computational modes will 
be exponentially dampened out. 

4. A few remarks 

• The proposed variants have improved numerical features in both numerical accuracy and computational effi- 
ciency compared to the original MCV3 scheme. 
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• In the new schemes, only the continuity of flux function is required at the cell boundaries where the constraint on 
the derivative of flux is not required anymore. This makes the schemes directly applicable to any quadrilateral 
and xahedral mesh. 

• Schemes with more solution points and higher order accuracy can be devised by the same spirit. 
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